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Abstract 

Through  a series  of  cases,  we  investigate  the  possibilities  of  the  shallow-water 
solver  Delft3D-Flow  of  simulating  the  evolution  of  (quasi-)2D  turbulence  in 
shallow  water  subjected  to  internal  and  external  friction  and  forcing. 

This  paper  presents  the  simplest  case,  namely  the  2D-DNS  of  laboratory 
experiments  of  freely-evolving  2D  turbulence,  initiated  by  a rake  in  a shallow 
fresh-water  layer  on  top  of  a salt-water  layer  in  a square  1*1  m2  basin  (Maassen 
2000).  Tabeling  et  al.  (1991)  report  similar  experiments  in  a fluid  with  a free 
surface  but  initiated  by  vortices  counter-rotating  in  a chessboard  arrangement. 
Likewise,  our  depth-averaged  free-surface  simulations  are  initiated  by  random  as 
well  as  by  chessboard  vortices.  We  compare  the  temporal  evolution  of  the 
simulated  vorticity  field  in  a viscous  fluid  with  observations  as  well  as  with 
simulations  of  Clercx  et  al.  (1999)  dedicated  to  incompressible  2D  turbulence 
simulation  with  a rigid  lid.  Theirs  and  our  simulations  agree  with  the 
experimentally  observed  evolution  of  vorticity  into  just  two  vortices  with 
opposite  rotation.  All  simulations  neglect  the  friction  at  the  density  interface  and 
exhibit  lesser  decay  of  kinetic  energy  than  observed  in  the  experiments  of 
Maassen  (2000). 

1.  Introduction 

Based  on  the  hydrostatic-pressure  assumption,  the  shallow-water  solver  Delft3D- 
Flow  has  been  extensively  used  for  modeling  the  depth-averaged  or  3D  flows  in 
civil-engineering  applications. 

For  better  assessment  of  ship  traffic,  structural  stability,  sediment  transport, 
dredging  operations  and  algae  blooms,  there  appears  to  be  a growing  interest  in 
simulating  more  details  of  the  flow  such  as  the  temporal  and  spatial  pdpf’s  of 
horizontal  velocity,  bed-shear  stress,  mixtures  of  dissolved  or  suspended 
constituents  etc. 

Rather  than  developing  2D  turbulence  closures,  we  prefer  resolving  and 
simulating  most  of  the  horizontal  flow  patterns.  Although,  practical  feasibility 
demands  that  we  maintain  modeling  the  3D  turbulence  in  the  orthodox  manner 
using  partial-slip,  bed  friction  and  3D  turbulence  closures  for  boundary-layer 
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type  of  flows.  We  define  the  latter  approach  as  Horizontal  Large  Eddy 
Simulation  (HLES)  acknowledging  the  particular  properties  of  shallow-water 
flows.  Recently,  this  HLES  approach  has  been  validated  against  observations  of  a 
shallow-water  mixing  layer  (Kernkamp  & Uittenbogaard,  2001). 

The  essential  question  is  whether  the  available  general-purpose  shallow-water 
solver  is  suitable  for  HLES,  particularly,  its  numerical  dissipation  and  the 
accuracy  constraints  for  simulating  2D-turbulence  in  free-surface  flows.  This 
paper  deals  with  these  fundamental  questions  by  simulating  laboratory 
experiments  of  2D  turbulence  freely-evolving  in  a square  basin. 

2.  The  shallow-water  solver  Delft3D-Flow 

In  depth-averaged  mode,  Delft3D-Flow  solves  the  following  depth-averaged 
shallow-water  equations  (SWE)  for  mass  conservation  and  for  the  depth- 
averaged  horizontal  velocity: 

+ V • (hu)  — Oand^=-  + u ■ Vw  + gV £ = vV2 u + 7^- h~lc f\u\u.  (1) 

S-)  t W ^ J I I 


In  (1),  u=(u,v)  is  the  depth-averaged  horizontal  velocity  vector,  h=d+£  the  total 
depth,  d the  still-water  depth  and  £ the  free-surface  position,  both  referring  to 
some  horizontal  reference  plane,  cf  the  bed-friction  coefficient,  v the  kinematic 
viscosity  and  g the  magnitude  of  gravitational  acceleration.  The  bed  friction  term 
is  the  depth  average  of  the  force  by  the  Reynolds-shear  stress  of  3D  turbulence 
without  wind.  The  force  vector  T is  due  to  subgrid-scale  stresses.  Bed  friction 
and  T are  omitted  in  this  paper  but  applied  by  Kernkamp  & Uittenbogaard 
(2001)  for  simulating  a shallow- water  mixing  layer.  The  term  gV  £ represents 
the  horizontal  gradient  of  pressure  based  on  the  hydrostatic-pressure  assumption. 
On  a staggered  (Arakawa-C)  grid,  the  horizontal  terms  in  (1)  are  formulated  in 
curvilinear  orthogonal  co-ordinates  but  in  this  paper  we  consider  just  Cartesian 
(x,y)-(u,v)  grids. 

The  SWE  (1)  are  time  integrated  by  ADI  of  which  the  first  half  time  step,  from 
t=nAt  to  (n+‘/2)At,  begins  by  solving  the  discretised  v-momentum  equation: 

+ unD?upw  t"+* }+  vnD2"upw  i"+4+  g£>f  cn'ri(r}= 

\At  x 1 J y 1 J y b J (2) 

2 vDf cmrl  {p  f cmrl  {v " }}+  2v  D ^ cn,rl  \p cmrl  {v " }} 


where  Da  is  a first-order  difference  operator  with  respect  to  a horizontal 
Cartesian  co-ordinate  a=x, y.  For  the  advection  terms  in  (2),  Da  is  2nd  order 
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upwind  but  2nd  order  central  for  pressure  as  well  as  stresses.  Eq.  (2)  is  implicit 
for  vn+,/\  solved  by  red-black  Jacobi  while  pressure  and  twice  the  Laplacian 
vV2v  are  integrated  explicitly.  Instead,  the  discretised  u-momentum  (3)  reads: 


n + \,p  + \ r 

ll  - - u 

Tai 


«+f.P+ 1 r.2" 

• + u - D, 


cntrl 


{un}+vn+b2;™r,i 


1=0. 

1 J I J 

n 2 


(3) 


Using  (3),  the  new  iterate  hn+-'p~\(n+2'p~]  t with  iteration  level  p,  is  expressed 
into  old-time  level  variables  as  well  as  into  £’”+I'-p+l  and  subsequently  substituted 
into  the  following  discretised  conservation  equation: 


r n+-j./’+l  _ ^ n 


4- A t 


+ D 


J.  M4 
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2"d cntrl  n+Lp+1 


l—./J  + l 


'+Dfcmrl{hnvn}=  0, 


(4) 


with  h centrally-averaged  to  (u,v)  points  yielding  a tri-diagonal  system  in  x- 
direction  that  is  solved  by  the  Thomas  algorithm  and  subsequently  iterated  on 

(p)  to  convergence.  Note  that  the  ratio  hn"'1'p  / hn+1'p+l  in  the  last  term  of  (3) 
ensures  mass-conservation  at  every  p-iteration.  In  the  subsequent  half-time  step, 
the  u-momentum  equation  is  formulated  and  solved  as  (2)  and  then  the  inviscid 
momentum  (3)  for  v is  coupled  to  (4)  while  (3)  is  formulated  implicitly  in  v- 
component  and  in  y-direction. 

For  the  full  time  step,  the  coupled  advection-conservation  schemes  are  second- 
order  in  time  with  fourth-order  dissipation  in  space  and  unconditionally  stable 
(Stelling,  1984).  The  viscous  force  is  integrated  explicitly,  invoking  a mild  time 
step  limitation. 

For  DNS  or  LES  applications,  a disputable  disadvantage  is  that  the  cyclic 
combination  of  2nd-order  upwind  and  2nd-order  central  advection  in  (2-4)  are 
neither  strictly  momentum  nor  strictly  energy  conserving.  On  the  other  hand, 
however,  the  tendency  of  creating  velocity  wiggles  in  inviscid  simulations  using 
the  energy-conserving  scheme  of  Arakawa  & Lamb  (1981)  appear  to  invoke 
more  dissipation  in  simulations  when  viscous  stresses  are  included  (Holm, 
1995).  Further,  the  semi-implicit  coupling  of  the  pressure  g£  to  the  momentum 
equations  (2-3)  as  well  as  to  the  conservation  equation  (4)  tends  to  conserve  the 
energy  stored  in  the  free-surface  or  compressible  velocity  modes  induced  by  the 
velocity-advection  operator.  Depending  on  the  flow  and  on  the  advection 
scheme,  the  latter  may  convert  incompressible  energy  modes  into  compressible 
energy  modes  that  are  subsequently  removed/dissipated  by  pressure 
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correction/projection  schemes  that  are  customary  in  fixed-volume 
incompressible  DNS/LES  solvers. 

This  section  is  concluded  by  introducing  the  time-step  limitations  based  on 
accuracy,  rather  than  stability,  of  advection  and  of  the  long-wave  or  barotropic 
(BT)  velocity  modes  on  a square  mesh  with  grid  size  Ax: 


&a=°a- 


A*  . _ <?BT  Ax  At, 

IT"  ‘272  7^7 


■ = 2V2  — 


■Jgh 


max 


BT 


<7bt  max 


0) 


(5) 


The  first  criterion  originates  from  the  semi-implicit  x-advection  term  in  (3)  and 
it  demands  aA<4  (Stelling,  1984,  p.  165).  Benque  et  al.,  (1982)  proposed  the 

second  criterion  with  aBx<4\/2  for  avoiding  aliasing  in  wave  propagation  using 
ADI  and  staggered  grids  e.g.  over  bed  topography.  For  the  temporal 
representation  of  advected  and  deforming  vortices  we  believe  aA<O(0.3)  would 
be  a proper  space-time  consistency  criterion.  For  Froude  numbers 
U I g4h  > 0(0.15)  , using  the  limitation  o BT  < 4^2  , the  latter  criterion 

coincides,  i.e.  AtA=AtHT.  However,  smaller  Froude  numbers  are  typical  for  our 
applications  and  then  the  accuracy  of  simulating  BT  modes  limits  the  time  step. 

3.  Some  considerations  of  2D  turbulence  with  a free  surface 

For  an  overview  of  properties  of  2D  turbulence,  such  as  the  inverse  energy 
cascade  related  to  vortex  merging,  we  refer  to  e.g.  (Lesieur,  1997).  For  an  closed 
basin  with  water  volume  V,  this  section  considers  briefly  the  possible  coupling 
of  2D  turbulence  (vortical)  motions  to  free-surface  (BT)  motions.  The 
frictionless  SWE  (1)  then  yields  the  following  conservation  property  (Arakawa 
& Lamb,  1981): 

j^\\h{\uu  + \gh-gd}dxdy  = 0,  (6) 


where  the  first  term  represents  the  kinetic  energy  (KE)  and  the  sum  of  second 
and  third  term  the  potential  energy  (PE).  For  a closed  basin  with  horizontal  bed 
(constant  d),  the  third  term  cancels  due  to  mass  conservation.  In  addition,  from 
the  frictionless  SWE  (1)  also  follows 


d2C 

Dt2 


ghV2C  = h{DijDij 


DC  <? 

— + u-v, 

Dt  dt 


(7) 
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with  horizontal-strain  rate  tensor  Du  based  on  u and  vertical  vorticity  component 
coz.  From  the  frictionless  SWE  also  follows  the  conservation  of  potential 
vorticity  (PV)  / h: 


D fo)^ 
Dt\h  , 


= 0. 


(8) 


A priori  (6-8)  suggest  that  2D  turbulence  in  free-surface  flows  differs  from  2D 
turbulence  in  rigid  fluid  volumes.  For  example,  the  LHS  of  (7)  represents  the 
long-wave  propagation  of  a surface  perturbation  £ associated  with  PE  in  (6) 
whereas  (6)  shows  that  PE  is  reversibly  exchanged  with  KE.  Further,  the  RHS  of 

(7)  acts  as  source  of  surface  perturbations  C,  of  which  the  vertical  vorticity  obeys 

(8) .  Note  that  (7)  is  similar  to  Lighthill’s  equation  for  sound  generation  by 
turbulence.  The  RHS  of  (7)  equals  the  Weiss-function,  see  e.g.  (Basdevant  & 
Philipovitch,  1994),  where  the  last  term  is  due  to  the  vertical  strain  rate  in  long 
waves. 

We  conclude  that  a priori  deviations  may  exist  between  the  evolution  of  2D 
turbulence  simulated  with  a rigid  lid,  such  as  by  Clercx  et  al.  (1999),  or  with  the 
SWE  (1).  Nevertheless,  detailed  analysis  (Vossen,  2000)  of  our  SWE 
simulations  as  well  as  animations  of  vorticity  combined  with  surface  elevation 
show  that  (7)  can  be  approximated  well  by  the  rigid-lid  counterpart  with 
pressure  g^: 

V2  (g(T)  -Co\  - D(j  Djj . (9) 

4 Numerical  initialisation  of  2D  turbulence  in  a free-surface  fluid 

Clercx  et  al.  (1999)  solve  the  viscous  2D  stream  function  and  vorticity  equations 
using  288"  Chebyshev  polynomial  expansion  and  they  initiate  2D  turbulence  by 
random  642  Chebyshev  modes  that  in  general  do  not  obey  incompressibility.  The 
laboratory  experiments  of  Tabeling  et  al.  (1991)  start  with  2D  chessboard 
vortices.  In  both  cases,  our  simulations  initiated  with  a spatially  constant  £ show 
that  much  PE  is  generated  as  freely-propagating  waves,  see  LHS  of  (7).  For  a 
proper  comparison  with  experiments  we  must  initiate  the  velocity-surface  field 
more  carefully  as  follows.  Define,  at  t=0,  the  master  velocity  field  as  u0  either 
based  on  642  Chebyshev  modes  or  chessboard  vortices.  Next,  the  associated 
incompressible  field  u(x,0)  is  derived  from: 

u(xfi)  =w0  + VAwithV2)l  =-  V -u0. 


(10) 
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with  Vk  the  corrected  velocity  field  such  that  V • u (x,0)  = 0 holds  on  the 
computational  grid.  For  closed  boundaries  (10)  yields  less  KE  in  u(x,0)  than  in 
u0.  Likewise,  Clercx  et  al.  (1999)  report  a decrease  in  KE  at  their  first  time  step 
that  yields  V • w(x,  Af)  = 0 . With  u(x,0)  by  (10),  the  initial  surface  elevation 
£(x,0)  is  obtained  from  (9)  as  the  quasi-steady  approximation  to  (7).  Following 
this  initialisation,  our  SWE  simulations  then  evolve  gradually  in  time  and  space 
without  notable  wave-like  motions. 


5.  Investigation  and  estimation  of  numerical  dissipation 

In  view  of  (5),  the  appropriate  time  step  was  investigated  by  inviscid  simulations 
of  (1)  initialised  by  Clercx’ s random  642  Chebyshev  modes  and  the  procedure 
given  in  section  4.  Figure  1 presents  the  temporal  evolution  of  volume-integrated 
KE  in  a basin  of  width  W=lm  with  1 cm  still-water  depth  and  lu'l=4  mm/s  using 
1002  square  grid.  In  all  cases,  cta<0.24  holds  but  Figure  1 shows  that  only  if 

<7bt  < 4>/2  is  satisfied  then  near  conservation  of  KE  is  obtained.  Therefore, 

(JBT  < 4V2  is  applied  in  all  subsequent  simulations. 


For  practical  reasons  we  estimate  a numerical  viscosity  vnum,  equivalent  to  the 
kinematic  viscosity,  by 


v 

num 


1 


dE 


2®  dt 


withE  = -t-pjjh^u2  +v2^xdyand&  = ^p^ho)2,dxdy,  (11) 


although  the  advection  operators  yield  4th  order  dissipation.  Series  of  inviscid 
runs  with  variations  of  grids  NP*NP  (Np=  50,  100,  200,  500),  of  basin  size  (1- 
100m)  and  initial  fields  with  random  Chebyshev  polynomials  or  chessboard 
vortices  (1,  4,  16  and  100)  yield  the  following  approximation: 


- . N =JL 

8^  ’ p Ax 


(12) 


with  Np  proportional  to  the  resolved  band  width.  Despite  the  strong  temporal 
merging  of  vortices  (see  figure  5),  Figure  2 presents  an  example  of  the  marginal 
temporal  dependence  of  vnum  estimated  through  (11)  and  thus  suggests  its 
validity  as  estimator  for  numerical  dissipation. 

6.  Simulations  of  2D  turbulence  in  free-surface  experiments 
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The  upper  part  of  figure  5 presents  particle  tracks  observed  in  a 1*1  m2  reservoir 
(Maassen,  2000).  Although  the  top  layer  floats  on  a denser  salt-water  layer, 
Maassen  (2000)  notes  a significant  interfacial  friction  reducing  the  effective 
Reynolds  number  by  a factor  5.  The  lower  part  of  figure  5 presents  the  vorticity 
contours  simulated  by  the  shallow-water  solver  but  without  interfacial  friction. 
The  evolution  of  vorticity  patterns  is  similated  qualitatively  well  by  our  shallow- 
water  solver.  Due  to  the  uncertainty  in  modelling  the  interfacial  friction,  i.e.  Cf  in 
(1),  we  prefer  the  comparison  with  simulations  by  Clercx  et  al.  (1999)  who 
applied  a solver  dedicated  to  this  type  of  experiments. 

For  Reynolds  number  Re=2000,  Clercx  et  al.  (1999)  simulated  decaying  viscous 
2D  turbulence  with  no-slip  wall  conditions  but  without  bed  or  interfacial  friction. 
They  define  Re  and  the  temporal  scale  T by 


Re  = 


\W\u' 

v 


Uy+~ 


Re 

17AL 


(13) 


with  W the  width  of  the  square  basin.  For  resolving  viscous  boundary  layers 
with  a square  grid  the  last  expression  in  (13)  should  equal  unity  and  this 
expression  determines  that  Np=200  is  adequate  for  a 2D  DNS  with  our  shallow- 
water  solver.  Figure  3 includes  the  simulations  of  Clercx  et  al.  (1999)  solved  by 
2882  Chebyshev  modes  and  we  applied  their  initial  velocity  field  based  on  642 
Chebyshev  polynomials  but  corrected  to  an  incompressible  flow  through  (10). 
Figure  3 presents  the  relative  decay  of  volume-integrated  KE  as  well  as  volume- 
integrated  enstrophy,  as  defined  by  (11),  against  time  scaled  by  T defined  in 
(13).  Under  these  conditions,  the  overall  numerical  viscosity  (12)  of  the  shallow- 
water  solver  is  estimated  to  be  about  0.2v.  Nevertheless  the  decay  simulated 
with  the  shallow-water  solver  is  comparable  or  even  less  than  simulated  by 
Clercx  et  al.  (1999). 


7.  Conclusions  and  discussion 

We  conclude  that  the  general-purpose  shallow-water  solver  Delft3D-Flow  is 
capable  of  simulating  vortex  merging  and  the  decay  of  2D  turbulence  in  a 
viscous  fluid  without  notable  numerical  diffusion  (figure  3).  From  a simulation 
point-of-view,  the  most  notable  aspect  of  the  free  surface  in  2D  turbulence  reads 
as  follows.  Animations  of  vorticity  and  surface  elevations  show  that  most  of  the 
free-surface  elevations  are  induced  by  the  Weiss  function,  RHS  in  (9).  Despite 
that  the  free-surface  elevations  carry  marginal  available  potential  energy, 
compared  to  kinetic  energy,  their  evolution  must  be  carefully  simulated  else 
kinetic  energy  is  generated  artificially  (Figure  1).  The  latter  demands  for  the 
barotropic  Courant  number  aBT  / 4 42  < 1 (figure  1 ) and  it  is  very  restrictive 
compared  to  other  accuracy  and  stability  conditions  for  the  shallow-water  solver 
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Delft3D-Flow.  The  reason  for  the  constraint  on  this  typical  wave-propagation 
condition  is  not  clear.  We  speculate  that  the  Poisson  equation  (9)  is  essential  and 
needs  to  be  solved  accurately.  In  SWE,  however,  the  solution  of  (9)  is  obtained 
through  (7)  and,  if  cBt  is  too  large,  errors  in  the  rapid  propagation  of  surface 
perturbations  spoil  the  approximation  to  (9). 
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Figure  1.  Decay  of  kinetic  energy  1002  grid,  1*1*0.01  m3,  lu'l=4  mm/s,  initial  random  field,  for 
<JBT  / 4-\/2  = 9.4  (6) ; 4.8  (5) ; 3.9  (4) ; 1 .9  (3)  ; 0.94  (2) ; 0.47  (1). 
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Figure  2.  Numerical 
500x500  (4)  grid  cel 
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Figure  3.  Relative  decay  of  volume-integrated  kinetic  energy  and  enstrophy,  see  (11)  of  2D 
turbulence  initiated  by  a random  field  at  Rc=2000  and  T=125  s time  scale.  Comparison  between 
the  shallow-water  solver  (Delft3D)  using  2002  square  grid  at  ca=0.0I2  and  anTi 4^/7  =0.47,  sec 
(5),  with  simulations  in  (Clercx  et  al.,  1999). 


Figure  4.  From  left  to  right  t=  0 ; 15  ; 20  ; 30  minutes.  Evolution  of  4*4  chessboard  vortices  on  a 5002 
grid  of  100*100*lm’ ; lu'l=0.4  m/s  ; inviscid  shallow-water  simulation. 
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Figure  5.  Particle  tracks  in  the  free-surface  water  layer  on  top  of  a dense  salt  layer  10 
seconds  (up-left),  1 minute  (up-right),  5 minutes  (down-left)  and  50  minutes  after 
passage  of  a raKE  (Maassen,  2000).  Below,  iso-contours  of  vertical  vorticity,  simulated 
by  Delft3D  for  conditions  of  figure  3 and  similar  to  observations  above. 
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